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ABSTRACT 

We present the first results from three-dimensional spherical shell simulations of magnetic dynamo action 
realized by turbulent convection penetrating downward into a tachocline of rotational shear. This permits 
us to assess several dynamical elements believed to be crucial to the operation of the solar global dynamo, 
variously involving differential rotation resulting from convection, magnetic pumping, and amplification of 
fields by stretching within the tachocline. The simulations reveal that strong axisymmetric toroidal magnetic 
fields (about 3000 G in strength) are realized within the lower stable layer, unlike in the convection zone 
where fluctuating fields are predominant. The toroidal fields in the stable layer possess a striking persistent 
antisymmetric parity, with fields in the northern hemisphere largely of opposite polarity to those in the southern 
hemisphere. The associated mean poloidal magnetic fields there have a clear dipolar geometry, but we have not 
yet observed any distinctive reversals or latitudinal propagation. The presence of these deep magnetic fields 
appears to stabilize the sense of mean fields produced by vigorous dynamo action in the bulk of the convection 
zone. 

Subject headings: convection - MHD - Sun: magnetic fields - turbulence 



1 . ELEMENTS OF SOLAR DYNAMO 

The solar global dynamo responsible for the observed 22- 
year cycles of magnetic activity is likely to require several dy- 
namical processes operating at differing sites within the con- 
vection zone and below its base. The primary elements were 
first discussed within the conceptual "interface dynamo" pro- 
posed by Parker (1993) and then explored by Charbonneau 
& MacGregor (1997). The key dynamical processes involve 
(a) the generation of magnetic fields by the intense turbu- 
lence influenced by rotation within the deep solar convection 
zone, (b) the transport or pumping of these fields downward 
into the tachocline of shear at the base of this zone, (c) the 
stretching of fields there by differential rotation to form strong 
toroidal fields, and (d) the magnetic buoyancy instability of 
such structures leading to field loops ascending toward the 
surface. Fully self-consistent magnetohydrodynamic (MHD) 
simulations of the complete solar global dynamo are not yet 
feasible, given the vast range of dynamical scales involved 
in these turbulent processes. However, many of the elements 
have now begun to be studied singly with reasonable fidelity 
using three-dimensional numerical modelling, often turning 
to localized planar domains to obtain sufficient spatial resolu- 
tion within the simulations (see reviews by Ossendrij ver 2003, 
Fan 2004, Charbonneau 2005, Miesch 2005). 

Yet some aspects require using full spherical shell geome- 
tries in the modelling. Early simulations of dynamo processes 
in spherical shells provided insights into the coupling of con- 
vection, rotation and magnetism, but were limited by spatial 
resolution (e.g., Gilman 1983; Glatzmaier 1985a, b). Recent 
studies of the interaction of turbulent convection with rota- 
tion (e.g., Miesch et al. 2000; Brun & Toomre 2002; Mi- 
esch, Brun & Toomre 2006) now yield global differential ro- 
tation profiles in close accord with helioseismic deductions 
(e.g., Thompson et al. 2003) in the bulk of the convection 
zone. Similarly, related MHD modelling of dynamo processes 



in such deep shells of turbulent convection has revealed that 
strong magnetic fields are produced without diminishing the 
differential rotation significantly (Brun, Miesch & Toomre 
2004, hereafter BMT04). Fluctuating magnetic fields with 
strengths of order 5000 G are realized that possess complex 
structures, with radial fields concentrated in downflow lanes 
and toroidal fields appearing as twisted ribbons extended in 
longitude. However, the associated mean fields are relatively 
weak and do not exhibit the systematic latitudinal propaga- 
tion or periodic polarity reversals seen in the Sun. Whereas 
BMT04 only considered dynamical element (a), we report 
here on also incorporating elements (b) and (c) into global 
simulations of dynamo action achieved by turbulent convec- 
tion able to penetrate downward into a tachocline of rotational 
shear. Such modeling is the next step in refining our intuition 
about the operation of the interface dynamo. 

2. CONVECTIVE SHELL WITH PENETRATION AND FORCED 
TACHOCLINE 

We conduct three-dimensional nonlinear MHD simulations 
of convection and dynamo action in a rotating spherical shell 
using the anelastic spherical harmonic pseudospectral code 
ASH (see Clune et al. 1999 and BMT04). ASH is based upon 
a large-eddy simulation (LES) approach involving subgrid- 
scale (SGS) modelling of unresolved turbulent processes. Our 
studies here extend the modelling of BMT04 by allowing 
penetration of the convection downward into a region of sta- 
ble stratification, where the pronounced differential rotation 
maintained by the convection is forced to vanish, thereby es- 
tablishing a tachocline of shear. Our computational domain 
extends over 0.62-0.96 R, where R is the solar radius, span- 
ning the bulk of the solar convection zone and part of the ra- 
diative interior below. Solar values are used for the rotation 
rate and luminosity, and the initial stratification is obtained 
from a one-dimensional solar structure model. We adopt a 
softer subadiabatic stratification in the lower stable zone to 
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ease our numerical resolution of internal gravity waves pro- 
duced by the overshooting motions (see Miesch et al. 2000). 

The tachocline as revealed by helioseismology is a nar- 
row transition boundary layer between the differentially ro- 
tating convection zone (fast equator, slow poles) and the uni- 
formly rotating radiative interior. Its discovery motivated the 
interface dynamo paradigm. The detailed structure of the 
tachocline may be determined variously by anisotropic turbu- 
lent mixing processes, magnetic stresses, and production of 
internal gravity waves; the equilibration times likewise con- 
tinue to engender debate (see review by Miesch 2005). De- 
spite uncertainties, a tachocline of rotational shear is crucial 
to element (c), and thus we seek to impose one here in two 
complementary ways. First, we introduce a drag force upon 
the axisymmetric velocities (relative to our rotational frame) 
to force them to vanish at the base of the computational do- 
main. We accomplish this smoothly with a hyperbolic tan- 
gent (of width 0.01 R centered at 0.66 R in the stable layer) 
so that motions within the bulk of the convective envelope 
are unimpeded. Without this forcing, the differential rota- 
tion established self-consistently within the convection zone 
would imprint itself upon the radiative zone through viscous 
and thermal diffusion. Second, we impose within the over- 
shooting region a weak latitudinal entropy variation to em- 
ulate the coupling between the convective envelope and the 
radiative interior through thermal wind balance within the so- 
lar tachocline. This is achieved by thermal forcing, involving 
a monotonic increase in temperature of about 6 K from equa- 
tor to pole at the base of convection zone (over a width of 
0.02 R), much as discussed by Miesch, Brun & Toomre (2006, 
hereafter MBT06). The thermal forcing aids in achieving a 
strong latitudinal differential rotation within the convection 
zone, which otherwise would be diminished by coupling to 
the radiative interior. This coupling is much stronger in our 
simulations than in the Sun because the viscous and thermal 
diffusivities are far larger and the overshoot region is wider. 

The domain boundaries are assumed to be impenetrable and 
free of viscous stresses. The lower boundary is assumed to be 
a perfect conductor and the radial entropy gradient is fixed 
there. At the upper boundary the magnetic field within the 
domain is matched to an external potential field and the en- 
tropy is fixed. The simulations were initiated by adding a 
weak seed field to a progenitor penetrative hydrodynamic sim- 
ulation. The seed field was toroidal and confined to the bulk 
of the convection zone. For the present model, we employ a 
horizontal resolution of Ng = 512, N<j, = 1024, and a stacked 
Chebyshev expansion in the radial dimension with N r = 98. 
The effective SGS viscosity v is equal to 6 x 10 12 cm 2 s _I at 
the top of the domain and decreases with depth as p -1 / 2 , where 
p is the mean density varying by a factor of 46 across the 
domain. The thermal (n) and magnetic (rj) SGS diffusivities 
have a similar profile, with fixed Prandtl and magnetic Prandtl 
numbers of Pr = v/k = 0.25 and Pm = v jr\= 8 respectively. 
The large Pm here reflects unresolved turbulent mixing pro- 
cesses, and such a value allows efficient dynamo action with 
tractable numerical resolution. In the Sun, Pm based on mi- 
croscopic processes is much smaller, but the relationship be- 
tween it and the effective turbulent diffusivities is uncertain. 
At small Pm, the critical magnetic Reynolds number needed 
for dynamo action increases considerably (Boldyrev & Catta- 
neo 2004; Schekochihin et al. 2005), rendering simulations 
much more computationally demanding. The field strengths 
reported here are likely sensitive to the Pm chosen. 




FIG. 1 . — (a) Global view (in Mollweide projection) at one instant of radial 
velocity v r on a spherical surface at mid-depth (r = 0.84S) in the convection 
zone, with upflows bright and downflows dark. Equator is shown dashed. An- 
gular velocity Q averaged in time and in longitude, showing (b) contours in 
radius and latitude and (c) variation with proportional radius along specified 
latitudinal cuts. 

Strong differential rotation is established by the convection. 
Figure lb shows the mean angular velocity ft (relative to the 
rotating frame) as a contour plot in radius and latitude, and 
Figure lc its variation with radius along selected latitudinal 
cuts. Within the bulk of the convection zone, there is a promi- 
nent decrease in ft from equator to pole. Near the upper 
boundary, the angular velocity contrast Aft (between equa- 
tor and 60°) is 55 nHz or 13% of the frame rotation rate. In 
comparison to MBT06 that had no penetrative region, the Aft 
here is reduced and 57 is more constant on cylinders aligned 
with the rotation axis. Below the base of the convection zone 
(at 0.73/?) where the stratification is subadiabatic, Q adjusts 
rapidly to uniform rotation within a pronounced tachocline 
of shear. The overall fast equator and slow pole behavior is 
in keeping with the £1 profiles deduced from helioseismology 
(e.g., Thompson et al. 2003), but this model possesses more 
radial shear within the bulk of the convection zone than in the 
Sun. We have clearly realized a tachocline of strong radial and 
latitudinal shear, which plays a major role in the organization 
of large-scale magnetic fields. 

Prominent magnetic dynamo action is achieved within this 
simulation. Much as in BMT04, the seed fields are rapidly 
amplified, with strong fluctuating magnetic fields (~ 3000 G 
rms) and weaker mean toroidal and poloidal fields (of time- 
averaged strengths < 300 G) achieved within the bulk of the 
convection zone. Over the 2800 simulated days (about 100 
rotation periods) studied after the field strengths have largely 
equilibrated, the dynamo action is persistent, with overall 
magnetic energy sustained at about 40% of the convective ki- 
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FIG. 2. — Toroidal magnetic fields realized in the convective envelope and underlying radiative region. Mollweide projections at one instant of longitudinal 
field B,p on a spherical surface (a) at mid-depth in the convection zone (r = 0.84J?) and (b) in the stable zone (r = 0.67K). (c) Contours in radius and latitude of 
averaged in time (over an interval of 220 days) and in longitude. Field strengths are in G. 



netic energy. Figure 2a shows a snapshot of the longitudinal 
field Z?^ at mid-depth (r = 0.847?), revealing complex structure 
on many scales, with no evident polarity preferences. The ra- 
dial and latitudinal fields B, and Bg also possess intricate and 
highly variable structures, tracing the convective flow realized 
in the bulk of the convection zone. In the underlying stable re- 
gion (at r = 0.67/?), Figure 2b shows that the magnetic field 
has been decidedly organized by the rotational shear, with 
the longitudinal field there stretched into large toroidal struc- 
tures that extend around much of the domain. Equally striking 
is that Bff, fields show antisymmetric parity in the stable re- 
gion, with opposite signs to those in the northern and southern 
hemispheres. The organized nature of strong toroidal fields in 
the stable region is confirmed in Figure 2c showing the time- 
averaged axisymmetric B$. Here the opposite polarities of the 
mean toroidal fields are evident in the two hemispheres, as 
contrasted to the weak and patchy mixed-polarity structures 
within the convection zone. The time-averaged axisymmet- 
ric fields in the stable region attain strengths of order 3000 
G, or about ten times stronger than the mean B^ in the con- 
vection zone. Whereas fluctuating (non-axisymmetric) fields 
dominate in the convection zone, the magnetic energy in the 
mean toroidal field is about three times larger than the fluc- 
tuating magnetic energy within the tachocline. The strong 
toroidal field established in the stable region is accompanied 
by a largely dipolar poloidal field. 

Neither the organization of the magnetic field below the 
convection zone into predominantly axisymmetric toroidal 
fields, nor the strong parity selection exhibited here, appear 
to be transient effects. These attributes arose rather quickly 
(~ 200 days) after the introduction of the forced tachocline, 
and have persisted for as long as we have continued the cal- 
culations. A sense of the evolution and spatial distribution 
of the mean fields in our simulation is provided by Figure 3 
showing the radial and temporal variations of the axisymmet- 
ric toroidal field B^ sampled at latitude -30°. The strongest 
toroidal fields there are clearly attained in the stable region 
(r/R < 0.73), with some modulations in amplitude but no 
changes in polarity. Within the convection zone, both the 
strength and polarity of the weaker toroidal field are more 
variable. Further investigation reveals that the axisymmetric 
dipole field both in the convection zone and in the stable re- 
gion has not changed its overall polarity during the course of 
our simulation. This is in sharp contrast to the evolution of 



the mean dipole component in BMT04, which considered the 
convection zone in isolation, where the dipole flipped at ir- 
regular intervals of less than 600 days. This suggests that the 
presence of a reservoir of strong toroidal fields with persistent 
polarity in the stable region is serving to stabilize the mean 
poloidal field realized in the convection zone. 

4. BUILDING STRONG TOROIDAL FIELDS 

The predominantly axisymmetric nature of the magnetic 
fields in the stable region may be understood in two comple- 
mentary ways. First, the lack of significant non-axisymmetric 
motions there precludes the generation of strong fluctuating 
fields like those realized in the convection zone. Thus non- 
axisymmetric (fluctuating) fields must either diffuse in from 
the convection zone or be transported downward by over- 
shooting motions. Second, Spruit (1999) has argued that any 
such fluctuating fields will, in the presence of well-defined 
angular velocity gradients, be quickly expelled from the sys- 
tem through reconnection between neighboring magnetic sur- 
faces. He estimates that this process acts on a time scale 
rn ~ {3r 2 ir 2 /ri&q 1 ) 1 / 3 , with q a measure of the rotational 
shear. In the Sun, this estimate yields time scales of order 100 
years, but in our far more diffusive simulation, tq is about a 
year. Thus fluctuating fields pervading the radiative region are 
fairly quickly erased; they are not replenished from above on 
a comparable time scale. 

The strength of the mean toroidal field below the convection 
zone is likely strongly influenced by the enhanced diffusive 
terms and spatially extended tachocline in our simulation. In 
the kinematic regime, the strength of toroidal fields B, gener- 
ated in time t by stretching of poloidal fields B p due to radial 
shear is of order B, /B p ~ Av^f / Ar, with Av^ the radial dif- 
ferential rotation velocity across a length Ar below the con- 
vection zone. Thus the strength of the toroidal field depends 
both on the angular velocity contrast and on how abrupt that 
contrast is; a narrower tachocline might then yield stronger 
fields. This estimate suggests that in our simulation, a 10 G 
poloidal field could be stretched to yield 1 kG toroidal fields 
in about a year. The magnetic diffusion time scale at the same 
depth is about 3 years, so fields can be amplified before they 
are diffused away. The far smaller magnetic diffusion and 
possibly narrower tachocline in the real solar interior may 
well lead to even greater amplification of the toroidal field 
there. 




FIG. 3. — Temporal evolution and spatial distribution of axisymmetric 
toroidal magnetic fields. Variation with radius and time of at a latitude 
of —30° , with tall peaks (bright tones) corresponding to large positive am- 
plitudes. The strongest fields are realized in the stable region (r/R < 0.73), 
where the polarity of the fields is also remarkably stable over the 2500 days 
sampled here. 



The realization of antisymmetric toroidal field parity within 
the stable zone is a striking property. It appears to be a robust 
feature but its origins are currently unclear. The generation of 
such parity would be expected if a dipolar poloidal field were 
to be sheared by differential rotation, but no such dipole seeds 
were introduced. Rather, the initial seed field was toroidal 
and included both symmetric and antisymmetric components. 
The emergence of a dipolar mean field may be a consequence 



of large-scale self-organization processes such as those asso- 
ciated with the inverse cascade of magnetic helicity in MHD 
turbulence or the a-effect of mean-field dynamo theory. Such 
processes may occur within the convection zone where Cori- 
olis forces induce opposite kinetic helicity in the northern and 
southern hemispheres (as in BMT04). Poloidal fields thus 
generated may then be pumped into the tachocline where they 
are stretched into toroidal structures by the shear and ampli- 
fied. The reservoir of toroidal field thus established within the 
tachocline may then feed back on field generation in the con- 
vection zone, providing the seed for new poloidal field. We 
plan to assess each of these processes in future work. 

In summary, these are the first global three-dimensional 
MHD simulations of turbulent solar convection to allow pen- 
etration into a forced tachocline of rotational shear. We have 
thus made contact with several of the dynamical processes 
thought to be essential in the operation of the global solar dy- 
namo. The combined action of convection and an organized 
tachocline has here yielded magnetic fields that possess sev- 
eral distinctive properties. The significant toroidal fields re- 
alized beneath the convection zone, the antisymmetric parity 
displayed by those fields, and the persistence of a single po- 
larity for multiple years, are reminiscent of the highly orga- 
nized magnetism that appears as sunspots at the solar surface. 
Further work to test the robustness of these features has been 
initiated and will be reported. 
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